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Abstract 

We introduce an approach to exploit the existence of 
multiple levels of description of a physical system to 
radically accelerate the determination of thermody- 
namic quantities. We first give a proof of principle of 
the method using two empirical interatomic potential 
functions. We then apply the technique to feed infor- 
mation from an interatomic potential into otherwise 
inaccessible quantum mechanical tight-binding calcu- 
lations of the reconstruction of partial dislocations in 
silicon at finite temperature. With this approach, 
comprehensive ab initio studies at finite temperature 
will now be possible. 
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Our understanding of the detailed processes under- 
lying the behavior of condensed matter systems has 
progressed tremendously over the last fifteen years. 
Electronic structure techniques give insight into the 
electronic origins of phenomena, but most readily for 
fixed lattices at zero temperature. Modern simple, 
but realistic, interatomic potentials reveal the lattice 
processes underlying complex phenomena at finite 
temperature. One of the greatest challenges in con- 
densed matter theory is to understand the electronic 
origins of solid state phenomena at finite tempera- 
ture. Comparatively little progress has been made 
to date in this area because of the tremendous com- 
plexity of combining the large number of degrees of 
freedom needed to describe the electrons with the 
large number of ionic positions needed to describe the 
occupied regions of the lattice phase space at finite 
temperature. In this letter we show that there is an 
exact separation of the finite-temperature electronic- 
structure problem into two far simpler parts, a finite- 
temperature lattice part, which may be studied us- 
ing approximate coarse-grained lattice models, and 



the electronic structure part, which may be studied 
using highly accurate, fine-grained calculations. 

Empirical interatomic potentials do not treat the 
electronic degrees of freedom explicitly, but have 
proved successful in modeling the general behavior of 
materials ranging from insulators ]]| through semi- 
conductors ^ |j| and metals These models cap- 
ture basic interatomic behavior and their simplicity 
makes them well suited for the evaluation of thermal 
averages. However, because empirical models coarse- 
grain over the electronic degrees of freedom, they gen- 
erally do not provide sufficiently accurate information 
for quantitative predictive studies. 

Tight-binding models represent a next step in de- 
tail of description and reliability. These models in- 
clude the electrons explicitly, but restrict their wave 
functions to linear combinations of atomic orbitals 
H ||, 0. These potentials are therefore applicable 
over a wider range of phase space than interatomic 
potentials and provide certain electronic information. 
Tight-binding calculations are far more demanding 
than their empirical interatomic counterparts. Di- 
rect Monte Carlo or molecular dynamics calculations 
of systems at finite temperature are sufficiently de- 
manding that the development of specialized approx- 
imate techniques is an active area of research || . 

Ab initio calculations attempt to describe all rele- 
vant electronic degrees of freedom. They give a good 
description of the physics of condensed matter sys- 
tems over a wide range of phase space jo). At present, 
density functional based calculations represent the 
greatest level of detail at which extended crystalline 
defects may be studied. Because the calculation of 
thermal averages with Monte Carlo or molecular dy- 
namics methods requires the evaluation of many con- 
figurations to sample phase space fully, ab initio cal- 
culations of thermal properties such as free energy 
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differences 10 are tremendously demanding and in- 
frequently attempted. 

Finite temperature studies require a strategy for 
obtaining a proper ensemble. In general, the popu- 
lated regions of phase space consist of a very narrow 
surface. To explore this surface with uniform sam- 
pling is impossible for a complex system because the 
surface occupies such a miniscule volume of phase 
space. As an alternative, one may start from a point 
on the occupied surface and then make small steps 
to explore it. One then always explores relatively 
relevant points in phase space, but covering the en- 
tire surface then requires an enormous number of 
steps. This can be manifested as correlation between 
large numbers of consecutive samples or as becom- 
ing trapped in local energy minima. This small-step 
approach is the essence of the two standard methods 
for evaluation of thermal averages: molecular dynam- 
ics and Monte Carlo calculations based on the 
Metropolis algorithm ]l2|| . 

As an alternative, we propose to separate the prob- 
lem into two parts. First, we fully explore the relevant 
regions of phase space using a coarse-grained Hamil- 
tonian, such as an interatomic potential for which 
extensive calculations based upon one of the tradi- 
tional small-step approaches are tractable. Once the 
relevant regions are identified, we evaluate physical 
observables within the target fine-grained Hamilto- 
nian, which may include electronic information, and 
combine the results to obtain proper averages over 
the fine-grained ensemble. Below we demonstrate 
that modern interatomic potentials are sufficiently 
close to their ab initio counterparts that thermal 
averages may be computed to within the accuracy 
of density functional theory by performing ab ini- 
tio calculations on a very limited number of samples 
drawn from the interatomic ensemble. Our approach 
takes optimal advantage of the radical separation in 
computational time scales which interatomic poten- 
tials and quantum mechanical calculations exhibit. 
Rather than employing ad hoc non-physical Hamil- 
tonians in a spirit similar to umbrella sampling |l3| . 
we propose the use of physical, albeit coarse-grained, 
Hamiltonians. 

Corrected Ensemble Averages — The first phase of 
our approach is to run a large-scale exploration of 
phase space with a simplified model, using an appro- 
priate, traditional thermodynamic approach. Draw- 
ing samples from the resulting configurations allows 
the more demanding model to be applied to truly 
uncorrelated, properly distributed points in phase 
space. So long as the two Hamiltonians are suffi- 
ciently correlated, we are assured that the samples 



selected are physically relevant to ensemble averages 
over the fine-grained Hamiltonian. 

To correct for the differences in the thermal distri- 
butions of the two Hamiltonians, each sample i se- 
lected from the initial simulation must be weighted 
by its relative probability in the two ensembles, the 
corrective Boltzmann factor 
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Here, /3 is the inverse temperature, and HI and H[ 
are the energy of configuration i within the coarse- 
grained and fine-grained models, respectively. Once 
given the corrective Boltzmann factors, the average 
of any observable O in the ensemble of H' is 
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Free energies — One common approach to the cal- 
culation of free energies is thermodynamic integration 
|l4j| , which gives the free energy difference between 
two macrostates characterized by the parameter val- 
ues A = and A = 1 as AF = f Q (^j) A dA, where 
the average is computed over the Boltzmann distribu- 
tion at fixed A. In practice, the integral is computed 
numerically by sampling a finite number of values of 
A. When A is a generalized coordinate, AF is the 
integral of the thermal average of the corresponding 
generalized force. 

Proof of Principle — To demonstrate the sound- 
ness of this approach we use two potentials which are 
sufficiently simple to allow direct comparisons using 
brute- force techniques. We have chosen two models 
which differ from each other at least as much as typ- 
ical interatomic potentials differ from density func- 
tional calculations. This allows us to address the is- 
sue of whether or not empirical potentials yield dis- 
tributions sufficiently close to ab initio distributions 
that (||) may be evaluated reliably using a tractable 
number of samples. 

To play the role of the coarse-grained Hamiltonian 
H c , we have chosen an early version of an interatomic 
potential for silicon which is still under development. 
The Stillinger- Weber model Q plays the role of the 
fine-grained Hamiltonian H' . The standard devia- 
tion in energy differences Hf — H c over thermally 
distributed samples provides a useful measure of the 
inter-Hamiltonian discrepancy. We found that, for 
the system we study below, this measure is the same 
(to within ten percent) for our two interatomic po- 
tentials as it is for the Stillinger- Weber potential and 
density-functional calculations. Later we shall apply 
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Figure 1: Unreconstructed (A = 1) configuration 
of the 30° partial dislocation (glide set) in silicon. 
Atoms in black are core atoms which bond together 
in pairs upon reconstruction. 

this measure to determine the relative suitability of 
different interatomic potentials in our procedure. 

The physical system to which we apply our method 
is the 30° partial dislocation (glide set) in silicon. 
Figure [l] shows the core of this dislocation in the 
unreconstructed state. In nature, the core atoms 
indicated in black in the figure move together in a 
period-doubling reconstruction to form dimers which 
saturate all dangling bonds and thereby minimize 
the energy of the dislocation. To illustrate our ap- 
proach, we calculate the free energy of this recon- 
struction at T — 930K, where some indirect experi- 
mental information is available JIB], |l6) . We use a pe- 
riodic supercell of 96 atoms containing a dislocation- 
antidislocation core at a separation of 14A. 

Figure ^| shows the cumulative changes in free en- 
ergy for TL C and TL^ (dashed and solid lines, respec- 
tively) as a function of A as we drive the dislocation 
core from its reconstructed (A = 0) to its unrecon- 
structed (A = 1) state. The figure also shows the 
statistical uncertainties remaining after the evalua- 
tion of twenty-five and fifteen million samples for TL C 
and W , respectively. As the figure illustrates, this 
number of samples is necessary to produce a reliable 
estimate of the final free energy difference within the 
Metropolis algorithm. Note that the final free energy 
change under TL C is much lower than for 7t* and that 
the energy profile for H c exhibits extraneous local 
minima and maxima. The challenge to our multiscale 
method is to correct for the far lower free energy and 
spurious local minima and maxima of TL C while using 
the same samples which led to the distorted curve for 



H c . 

The dash-dotted line Figure || displays the results 
we obtain with the multiscale approach. To produce 
these results, we first made five separate runs each 
drawing only one thousand independent samples at 
large intervals from the direct simulation under H. c . 
We then evaluate the generalized force using TL f and 
compute the average force at each value of A for each 
run according to (0). Finally, we find the mean and 
standard deviation among the forces of the five runs 
at each A and integrate the results numerically to 
give the free energy curve and uncertainties that ap- 
pear in the figure. The fact that the newer curve has 
quite similar statistical uncertainties to the previous 
curves while computed from three orders of magni- 
tude fewer samples underscores the fact that the vast 
majority of samples in the brute- force Metropolis ap- 
proach only serve to generate the Boltzmann distribu- 
tion but do not contribute significantly independent 
statistical information to the averages. If a sampling 
scheme better than the Metropolis algorithm were 
employed, the comparison would be less favorable for 
our approach. However, we note that our approach 
is direct, efficient and as general as the availability 
of suitable atomistic potentials. Figure || also shows 
that the averaging procedure (||) succeeds in eradicat- 
ing the spurious maxima and minima of the coarse- 
grained Hamiltonian. The final free energy of recon- 
struction determined using the multiscale approach 
is 0.711 ± 0.019 eV, in good agreement with the re- 
sult of the brute-force method, 0.712 ± 0.010 eV. We 
thus see that the low free energy difference under Ti. c 
is indeed properly compensated, even when working 
with a radically reduced number of samples. 

To obtain a statistical uncertainty in the free en- 
ergy difference of 0.043 eV, which is within the accu- 
racy of the best density functionals Q , the preceding 
calculation could be done with a single run, instead 
of five. This would require the evaluation of a total 
of only one thousand samples, quite feasible for ah 
initio work. 

While it is true that the corrective factors fluctu- 
ate exponentially with the total energy discrepancy 
Tif — TL C , whose standard deviation scales as the 
square root of the number of atoms in the system, the 
corrected average (||) does not fluctuate among runs 
nearly as widely. This is because (j2j) is a weighted 
average of the observable O, and, therefore, the fluc- 
tuations of O place an absolute upper bound on the 
fluctuations of (||). The system size dependence of 
the fluctuations of (||) thus ultimately becomes the 
system size dependence of the fluctuations of O. If 
O is a local or intensive parameter, as we have here, 
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Figure 2: Free energy of the silicon 30° partial dis- 
location (glide set) as driven from its reconstructed 
(A = 0) to its unreconstructed (A = 1) state. 

these fluctuations do not increase with system size. 

For runs of fixed length, increases in fluctuations in 
the corrective factors result in fewer values of O con- 
tributing to the final average. It thus becomes more 
difficult to detect the correlation between the observ- 
able and the inter-Hamiltonian discrepancy. How- 
ever, local observables, such as our generalized force, 
are largely uncorrelated with the global energy dis- 
crepancy. When these quantities are completely un- 
correlated, (|) reduces to < O >/= (l/iV) Ei i- 
This direct average yields a free energy of reconstruc- 
tion of 0.682 ± 0.008 eV, in excellent agreement with 
the exact result. The correlation between the gener- 
alized force and the energy discrepancy thus corre- 
sponds to the remaining 0.02 to 0.04 eV in the free 
energy, which the corrective Boltzmann factors cap- 
ture rather well within our present run length and 
system size. Numerical experiments on model ran- 
dom variables reveal that, for a fixed run length, the 
effect of correlation between a local observable and 
the energy discrepancy degrades slowly with system 
size, approximately as the square root of the number 
of atoms. As the present scale of one hundred atoms 
approaches the maximum size feasible for current ah 
initio calculations with extensive exploration of phase 
space, we do not expect this degradation to present a 
significant problem for some time. Much can be done 
on the scale of one hundred atoms with our approach 
in its present formulation. 

Multiscale model of dislocation cores at finite tem- 
perature — As a first truly multiscale application of 
our approach, we shall now demonstrate the use of 
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Figure 3: Multiscale calculation of quantum mechan- 
ical free energy of reconstruction 

interatomic potential functions, which do not deal 
explicitly with electronic degrees of freedom, to gen- 
erate samples for quantum-based tight-binding cal- 
culations, which do. The traditional Monte Carlo 
approach would demand the evaluation of millions of 
tight-binding configurations, requiring years of com- 
puting time. 

For the present calculations, we use the tight- 
binding Hamiltonian of Sawada [p| with the modi- 
fications proposed by Kohyama |7|]. This model pro- 
vides a relatively good description of the bulk, dimer 
and surface energetics of silicon. A limited number 
(about fifty) of preliminary runs on our target sys- 
tem shows that the standard deviation in energy dis- 
crepancy between this tight-binding model and the 
two interatomic potentials are 0.40 ± 0.04 eV and 
0.57 ±0.07 eV, with the Stillinger- Weber yielding the 
lower value. Accordingly, we use the Stillinger- Weber 
potential to explore the phase space of our system. 

Figure |^ presents the quantum mechanical tight- 
binding free energy profile of reconstruction calcu- 
lated with the multiscale approach. Our results are 
based on the energy and forces of only twelve hun- 
dred samples drawn from the Stillinger- Weber poten- 
tial. With the traditional Monte Carlo approach, a 
run of such length would produce at most a few inde- 
pendent points in phase space, while our calculation 
represents an extended exploration of the canonical 
distribution. Our tight-binding result for the free en- 
ergy of reconstruction is 0.53 ± 0.015 eV per broken 
bond in the unreconstructed dislocation. 

This value falls within limits placed by available ex- 
perimental information and lends support to mecha- 
nisms for dislocation mobility suggested by atomistic 
studies |19]. In these atomistic studies, the energy re- 



4 



quired to break the reconstructed bonds in the dislo- 
cation core contributes significantly to the activation 
energy for dislocation mobility. The Stillinger- Weber 
value for the reconstruction energy, 0.81 eV, however, 
is too high in that it leads to a prediction [[l8| for the 
dislocation mobility which is multiple orders of mag- 
nitude lower than observed experimentally Jl5[ . Our 
value for the free energy includes explicitly both ther- 
mal fluctuations and bonding effects. The fact that it 
is significantly lower than the Stillinger- Weber value 
tends to increase the predicted dislocation mobility 
and thus lends support to the mechanisms proposed 
in [|l9|. Our lower free energy leads also to a signif- 
icantly higher prediction for the equilibrium density 
of dangling bonds in the dislocation core. The total 
signal strength in electron paramagnetic resonance 
experiments fill , whose precise origin is difficult to 
interpret and may involve other defects, places an ab- 
solute upper limit of a few percent on the density of 
dangling bonds associated with the dislocations. Our 
free energy value corresponds to a dangling bond den- 
sity of ~ 0.1%, well below the experimental bound. 

Conclusions — There is tremendous benefit in sep- 
arating thermal studies of a system with highly de- 
tailed models into two parts: an exploration of phase 
space using a simpler coarse-grained Hamiltonian and 
use of a more detailed Hamiltonian to study the be- 
havior of the system at a limited number of well- 
chosen points in phase space. We have seen that 
the approach works well for the determination of the 
free energies of defects from quantum-based calcula- 
tions. As system size increases, the determination of 
such local defects is relatively stable. Global changes, 
such as phase transitions, in which the observable of 
interest is correlated with the entire system may re- 
quire a different approach. A thermal, quantum me- 
chanical description of the 30° partial dislocation core 
gives a free energy of reconstruction which is signif- 
icantly lower than previous Stillinger- Weber values 
and, thus, leads to a more consistent view of dis- 
location mobility in silicon. Finally, the number of 
samples required to apply the multiscale approach to 
systems of approximately one hundred atoms is suf- 
ficiently low as to make its application to ab initio 
calculations attractive. 
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